Anderson Localization in Metamaterials with Compositional Disorder 
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We consider one-dimensional periodic-on-average bi-layered models with random perturbations in 
dielectric constants of both basic slabs composing the structure unit-cell. We show that when the 
thicknesses d a and dt of basic layers are essentially nonequal, d a ^ d^, the localization length Li oc 
is described by the universal expression for two cases: (a) both layers are made from right-handed 
materials (the RH-RH model), (b) the a layers are of a right-handed material while the b layers 
are of a left-handed material (the RH-LH model). For these models the derived expression for Li oc 
includes all possible correlations between two disorders. However, when d a = dt the RH-LH model 
exhibits a highly non-trivial properties originated from inhomogeneous distribution of the phase of 
propagating wave, even in the case of white- noise disorder. We analytically show that in this case the 
localization length diverges in the conventional second order in perturbation parameters. Therefore, 
recently numerically discovered anomalies in L; oc are due to the next order of approximation. On 
the other hand, for the RH-RH model the general expression for Li oc remains valid for d a = db as 
well. 

PACS numbers: 72.15.Rn, 42.25.Dd, 42.70.Qs, 

Keywords: Anderson localization, Photonic crystals, Metamaterials 



-a 
c 
o 

> 
^1- 



X 



I. INTRODUCTION 

Due to progress in nano- and material science, the 
study of wave propagation and electron transport in pe- 
riodic one-dimensional (ID) systems has attracted much 
attention (see, e.g., Ref. 1 and references therein). The 
systems of particular interest are bi-layer structures in 
optics and electromagnetics, or semiconductor superlat- 
tices and arrays of alternating quantum wells and barriers 
in electronics. The interest to this subject is due to a pos- 
sibility to create structures with prescribed properties of 
transmission and/or reflection. New perspectives in this 
direction are related to unconventional optic properties 
of metamaterials. 

One of the practical problems is the influence of dis- 
order that cannot be avoided in real experiments. The 
disorder can be originated from fluctuations of the thick- 
ness of layers (positional disorder) or from variations of 
the medium parameters, such as permittivity and perme- 
ability for electromagnetic waves or the barrier height for 
electrons (compositional disorder). Typically, the disor- 
der is treated as an unwanted effect, however, recently 
it was understood that it can be a promising candidate 
for targeted manipulation of transport properties. In- 
deed, the correlations in the disorder may result in un- 
usual features of transport. In particular, it was shown, 
both theoretically [2-8] and experimentally [9-11], that 
specific long-range correlations can significantly enhance 
or suppress the wave localization in a desired window of 
frequency 

As is well known, the transmission through any ID dis- 
ordered system is governed by Anderson localization (see, 
e.g., Refs. 8, 12, 13 and references therein). The principal 



concept of this phenomenon is that all transport charac- 
teristics depend solely on the ratio L/Li oc between the 
structure length L and localization length Li oc of eigen- 
states. Such a universal dependence manifests itself, for 
example, in the famous expression for the self-averaging 
logarithm of the transmittance Tl, 



(lnT L ) =-2L/L 



loc- 



(1.1) 



Here, the angular brackets (. . .) stand for averaging over 
the disorder. In agreement with the concept of single pa- 
rameter scaling, there are only two characteristic regimes 
in ID disordered structures, namely, the regimes of ballis- 
tic and localized transport. The ballistic transport occurs 
when the localization length L\ oc is much larger than the 
sample length L. In this case the samples are practically 
fully transparent since its average transmittance is close 
to unity, 



(T L ) a 1 - 2L/L 



loc 



for 



L <£L L 



loc- 



(1.2) 



In the opposite case when Li oc is much smaller than L, 
the ID disordered structures exhibit localized transport. 
In this case the average transmittance is exponentially 
small, 



(Tl) « cxp(-2L/L 



loc) 



for Li oc < L 



(1.3) 



As one can see, in the localization regime a ID disordered 
system perfectly (with an exponential accuracy) reflects 
quantum or classical waves because of strong localiza- 
tion of all eigenstates. This brief analysis shows how one 
can control the transport by manipulating the value of 
localization length in comparison with the system size. 

In spite of a remarkable progress in understanding the 
main features of the wave (electron) propagation through 
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random structures, the majority of studies are based on 
numerical methods [14-25]. Giving the important results 
the numerical approaches obviously suffer from the lack 
of generality being restricted by specific values of param- 
eters. As for the analytical results, they are mainly refer 
to simplest models with white-noise disorder [26-30] or 
to the correlated disorder with delta-like potentials in the 
Anderson [2] or Kronig- Penney models [4, 31]. 

In this paper we address the ID models with periodic- 
on-average bi-layered structures with weakly disordered 
parameters. In comparison with Refs. 32, 33 where the 
case of random thicknesses of slabs was considered, here 
we focus on the model with weakly disordered refractive 
indices of both slabs. The main attention is paid to the 
comparison between normal stacks (when both layers are 
of normal material) and mixed stacks (with alternating 
normal and meta- material layers). First, we review the 
recent study of the structures with different thicknesses 
of two basic layers [34, 35] . Here we derive the unique an- 
alytical expression for the localization length L/ oc , which 
is valid in a large range of model parameters, and can be 
applied to various physical realizations. The key point 
of our approach is that we explicitly take into account 
possible correlations within two disorders of each layer 
type as well as between them. Second, we extend our 
analytical analysis to the particular case of equal thick- 
nesses of basic layers, see Ref. 23. We show that in this 
specific case the standard perturbation theory fails if one 
of two basic stacks is made from a left-handed material. 
Specifically, the localization length diverges in the con- 
ventional second order in perturbation parameters, thus 
leading to anomalous properties of scattering. 



II. PROBLEM FORMULATION 

The model describes the propagation of electromag- 
netic waves of frequency to along an infinite array of two 
alternating a and b layers. Every kind of slabs is respec- 
tively specified by the dielectric permittivity e a ,b, mag- 
netic permeability fi a ,b, corresponding refractive index 
n a ,b = ^/£ a ,b^a,b, impedance Z a ,b = fJ>a,b/n a ,b and wave 
number k a b = Lon ab /c. We consieder two cases: (a) 
when both a and b slabs consist of a right-handed (RH) 
optic materials, and (b) when a slabs are of RH-material 
and 6 layers are of left-handed (LH) material. In what 
follows, the combination of RH-RH slabs is referred to 
the homogeneous stack whereas the array of RH-LH lay- 
ers is called mixed stack. 

As is known, for the RH medium all optic parameters 
are positive. On the contrary, for the LH material the 
permittivity, permeability and corresponding refractive 
index are negative, however the impedance remains pos- 
itive. Every alternating slab has the constant thickness 
d a or db, respectively, so that the size d of the unit (a, b) 
cell is also constant, d = d a + db- 

As was noted in Ref. 33, when the impedances of two 
basic a and b slabs are matched, Z a = Z b , the localization 



length diverges and the perfect transparency emerges, 
while a positional disorder itself persists. In the follow- 
ing, we analyze the effect of compositional disorder in a 
stack-structure whose unperturbed counterpart consists 
of perfectly matched basic a and b layers. Specifically, 
following Ref. 23, we admit that a disorder is incorpo- 
rated via the dielectric constants e a ,b only, so that 

£a(n) = [1 + r] a (n)} 2 , n a = 1, 

n a {n) = l + r] a (n), Z a {n) = [1 + r^n)]" 1 ; (2.1a) 

e b (n) = ±[1 + rjb(n)} 2 , /z& = ±1, 

n b {n) = ±[1 + Vb (n)}, Z b (n) = [1 + r 7fe (n)]- 1 .(2.1b) 

Here integer n enumerates the unit (a, b) cells. The upper 
sign stands for the RH material while the lower one is 
associated with LH medium. 

Without disorder, f] a . b (n) — 0, the homogeneous RH- 
RH structure is just the air without any stratification, 
while the mixed RH-LH array represents the so-called 
ideal mixed stack (e a = fj, a = 1, s b = Mb = —1) with 
perfectly matched slabs (Z a = Z b = 1). Therefore, both 
the unperturbed RH-RH and RH-LH stacks are equiva- 
lent to the homogeneous media with the refractive index 
n, thus resulting in no gaps in their linear spectrum, 

-/ - \d a ±d b \ 
K = u}n/c, n—— — . 2.2 

da+db 

The random sequences rj a {n) and r] b (n) describing 
the compositional disorder, are statistically homogeneous 
with zero average, {r] a ,b{n)) — 0, and binary correlation 
functions defined by 

(Va{n)va(n')) = o- 2 a K a (n - n') , (2.3a) 
(Vb(n)n b (n')) = a 2 K b (n - n') , (2.3b) 
(Va(n)vb(n')) = cr 2 ab K ab (n - n') . (2.3c) 

The averaging (...) is performed over the whole array or 
due to the ensemble averaging, that is assumed to be 
equivalent. The auto-correlators K a (r) and Kb(r) as well 
as the cross-correlator K a b(r) are normalized to unity, 
K a (0) = K b (0) = K ab (0) = 1, and decrease with an 
increase of the distance |r| = \n — n'\ between the cell 
indices n and n' . The variances a 2 and a 2 are obviously 
positive, however, the term a 2 b can be of arbitrary value 
(positive, negative or zero). We assume the composi- 
tional disorder to be weak, i.e. 

6 < 1, (ka,bd a ,b) 2 <ylb < 1; (2-4) 

that allows us to develop a proper perturbation theory. 
In this case all transport properties are entirely deter- 
mined by the randomness power spectra IC a (k), lCb(k), 
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and )C a b(k), denned by the standard relations 



il>' a b at the opposite boundaries of the same a or b layer, 



JC{k) = K(r)exp(-ikr), (2.5a) 

r— — oc 

K(r) = — f dklC{k)exp(ikr). (2.5b) 

271 " J-7T 



By definition (2.3), all the correlators K a (r), Kb(r) and 
Kabir) are real and even functions of the difference 
r = n — n' between cell indices. Therefore, the corre- 
sponding Fourier transforms JC a (k), JCb(k) and fC a b{k) are 
real and even functions of the dimensionless wave number 
k. It should be also stressed that according to rigorous 
mathematical theorem, the power spectrum JC(k) is non- 
negative for any real sequence. 

Within any of a or b layers the electric field of the 
wave i[)(x) exp(-iujt) obeys the ID Hclmholtz equation 
with two boundary conditions at the interfaces between 
neighboring slabs, 



(^ + kl b )^ a , b (x)=0, 



\dx 2 



(2.6a) 



ipa{xi) = i>b{xi), Ha Vafci) = H Vb(zi)- (2-6b) 



The x-axis is directed along the array of bi-layers per- 
pendicular to the stratification, with x — x, standing for 
the interface coordinate. 



III. BASIC EXPRESSIONS 

The general solution of Eq. (2.6a) within the nth (a, b) 
cell can be presented in the following form, 

i> a (x) = il) a {x an ) cos [k a (x - x an )\ 

+K 1 i ) 'a{ x an) SUl [k a (x - X an )] (3.1a) 

inside a n layer, where x an < x < x bn ; 

1p b (x) = 4>b (Xbn ) COS [k b (x -Xbn)] 

+k b 1 ip' b (x bn ) sin [k b (x - x bn )] (3.1b) 
inside b n layer, where x bn < x < x a (n+i) . 

The coordinates x an and x bn refer to the left-hand edges 
of successive a n and b n layers. Note that Xb n — x an = d a 
and x a ( n+1 ) — Xb n — d-b- The solution (3.1) gives a useful 
relation between the wave function ip a ^ and its derivative 



1pa(Xbn) = 1pa(Xan)cOS(fi a (n) 

+K 1 ip' a (x an )smp a (n), 

i>'a( X bn) = -k a 1pa(x a n)smlfi a (n) 

+ip' a {x an )cos^ a {n)\ (3.2a) 

A(x a (n+1)) = Ipb(Xbn) COS p b (n) 

+k b 1 ip' b (x hn ) sin^ & (n), 

i>'b{ x a{n+l)) = -kb4>b(Xbn) SUl (f b (n) 

+ip' b {xbn) cos <pb{n). (3.2b) 

The corresponding phase shifts ip a (n) and (pb( n ) depend 
on the cell index n via random refractive indices (2.1), 

!fia(n) = k a (n)d a = ip a [l + Va(n)], <p a = ud a /c; (3.3a) 

<p b (ri) = kb(n)d b = <pb[l + Vb(n)], <p b = ± wd(,/c. (3.3b) 

By combining Eqs. (3.2) with boundary conditions (2.6b) 
at Xi — x bn and Xi = x a ( n +i), one can write the recurrent 
relations for two opposite edges of the whole nth unit 
(a, b) cell, 

Qn+l = A n Q n + B n P n , Pn+l = —CnQn+D n P n . (3.4) 

Here the "coordinate" Q n and "momentum" P n are 

Qn = il>a{x a „) and P n = (c/u})i)' a (x an ), (3.5) 

with indices n and n + 1 standing for left and right edges 
of the nth cell. The factors A n , B ni C n , D n read 

A n = cos(p a cos(pb ~ Z^ZbsmipaSincpb, (3.6a) 

B n = Z a smip a coslp b + Z b cos(p a smlpb, (3.6b) 

C n = Z^ 1 sin <p a cos tpb + Z^ 1 cos <p a sin tp b , (3.6c) 

D n = cos [p a cos Lpb - Z a Z^ x sin (f a sin <p b - (3.6d) 

They depend on the cell number n due to random refrac- 
tive indices (2.1) that enter into both the impedances 
Z a ,b{n) and phase shifts <f a ,b{n). It is noteworthy that 
the recurrent relations (3.4) can be treated as the Hamil- 
tonian map of trajectories in the phase space (Q, P) with 
discrete time n for a linear oscillator subjected to time- 
depended parametric force. 

With vanishing disorder, n 0j b(n) = 0, the factors (3.6) 
do not depend on the cell index (time) n. Therefore, 
in line with map (3.4), the trajectory Q ni P n creates a 
circle in the phase space (Q,P) that is an image of the 
unperturbed motion, 



Qn+l = Qn COS 7 + P n Sm 7, 
Pn+l = -Qn Sin 7 + P n COS 7. 



(3.7) 



The unperturbed phase shift 7 over a single unit (a, b) 
cell is defined as 



7 = ipa + fb = u(d a ± d b )/c. 



(3.8) 
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This result is in a complete correspondence with the spec- 
trum (2.2) taking into account that the Bloch wave num- 
ber k — \l\/d. 

Having the circle (3.7), it is suitable to pass to polar 
coordinates R n and 9 n via the standard transformation 



Qn = Rn COS ( 



P n = R n sin ( 



(3.9) 



By direct substitution of Eq. (3.9) into map (3.7), one can 
see that for the unperturbed trajectory the radius R n is 
conserved, while the phase n changes by the Bloch phase 
7 in one step of time n, 



R 



n+l 



— Rn 



3 n+l 



7- 



(3.10) 



Evidently, the random perturbations, r] a , b (n) ^ 0, give 
rise to a distortion of the circle (3.10) that is described 
by the Hamiltonian map (3.4) with randomized factors 
(3.6). With the use of definition (3.9), one can readily 
rewrite this disordered map in the radius-angle presen- 
tation. The corresponding exact recurrent relations read 



p2 

fV+i 
Rl 



= {An + C 2 ) cos n + (B n + D n ) sin 0. 



(A n B n - C n D n )sm20 n , 

— C n + D n tan n 



tan 6. 



n+l 



A n + B n tan n 



(3.11a) 
(3.11b) 



Eqs. (3.11) constitute the complete set of equations in 
order to derive the localization length Li oc according to 
its definition via the Lyapunov exponent A [36, 37], 



Bloc 2 \ V Rn 



(3.12) 



Note that in Eq. (3.12) the averaging is performed along 
the trajectory specified by R n and n . 

Here we should emphasize the following. In the ideal 
mixed stack (e a = fi a = 1, s b = fi b = -1, Z a = Z b = 1) 
the wave spectrum (2.2) is singular. Specifically, when 
thicknesses d a and d b are equal, d a = d b , the phase 
velocity c/n diverges and the Bloch phase 7 vanishes. 
As a result, the circle (3.10) presenting the unperturbed 
map, degenerates into a point in the phase space (Q, P). 
Therefore, the perturbation theory has to be developed 
in a different way depending on the value, finite or van- 
ishing, of the Bloch phase 7. For this reason, in what 
follows we perform the separate analysis for d a — d b case 
when considering the RH-LH stack-structure. 



IV. BI-LAYER ARRAY WITH FINITE BLOCH 
PHASE 

In this section we assume the arbitrary relation be- 
tween slab thicknesses, d a and d b , for the homogeneous 
RH-RH bi-layer array, however, for mixed RH-LH stack- 
structure we assume only different thicknesses of basic 



slabs, d a ^ d b . In this case the Bloch phase has finite 
value, 7^0, and the localization length can be derived 
in the standard way already used in the previous studies 
[2-4, 31, 33, 34]. Specifically, we expand the coefficients 
(3.6) up to the second order in the perturbation param- 
eters r) a ^ b (n) <C 1. In doing so, one can expand only 
the factors containing the impedances Z a ^(n). As to the 
random phase shifts {p a ,b{n), they can be replaced with 
their unperturbed values <p a ,b, see definitions (3.3). This 
fact becomes clear if we take into account the conclusion 
of Ref. 33: The phase disorder contributes to the Lya- 
punov exponent only when the unperturbed impedances 
are different. The quite cumbersome calculations allow 
us to derive the perturbed map for the radius R n and 
angle n , 



n n+l 



I + Va(n)V a (0 n ) + Vb(n)V b (0 n ) 



R 2 

dL n 

+V 2 aHW a + if b {n)W b + Va(n)Vb(n)W ab , (4.1a) 

n +i -0 rl +7 = Va(n)U a (0 n ) + Vb (n) U b (0 rl ) . (4.1b) 

Here the functions standing at random variables r) a , b {n) 
are described by the expressions: 

V a (0 n ) = -2sm(p a am(20 n - <p a ), (4.2a) 

V b (0 n ) = -2sin<£ 6 sin(20„ -7- <^ ), (4.2b) 

W a = 2 sin 2 (p a , W b = 2 sin 2 tp b , (4.2c) 

W ab = 4 sin ip a sin ip b cos 7; (4.2d) 

U a {0 n ) = - sin ^ a cos(26»„ - ip a ), (4.2e) 

U b {0 n ) = - sin ip b cos(20„ - 7 - <p a ). (4.2f) 

In Eqs. (4.2) we keep only those terms that contribute to 
the localization length Bi oc in the first non-vanishing or- 
der of approximation. Note that in Eq. (4.1a) the factors 
V aib containing n are always multiplied by rj atb (n), there- 
fore, only linear terms in the perturbation are needed in 
the recurrent relation (4.1b) for the angle 0„. 

Now, in order to evaluate the Lyapunov exponents one 
has to substitute Eq. (4.1a) into Eq. (3.12) and expand 
the logarithm within the quadratic approximation in the 
perturbation parameters r} a ,b{n)- In this approximation 
one can treat the terms r] 2 (n), rj 2 (n) and r] a (n)r] b (n) as 
uncorrelated with factors (4.2) containing the angle vari- 
able n - It is important that performing the averaging we 
assume the distribution of phase n to be homogeneous, 
i.e., the corresponding distribution function p{6) = l/2n. 
One can show that this assumption is valid apart from 
the case 7 = 0, i.e. when we consider the ideal mixed 
RH-LH stack-structure with d a = d b . After some alge- 
bra we arrive at the final expression for the Lyapunov 
exponent, 

A = = \a 2 a lC a {2 1 ) sin 2 (p a + ^of/C 6 (27) sin 2 ip b 

+al b lCa b (2j) sin Lp a sin (p b cos (4.3) 
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As one can expect, the result (4.3) is symmetric with 
respect to the permutation of slab indices a <-> b. 

The expression (4.3) for the Lyapunov exponent A con- 
sists of three terms. The first two terms contain the auto- 
correlators K a {2 r y) or JCb(2j), they are originated from 
the correlations between solely a or solely b slabs, respec- 
tively. The third term depends on the cross- correlator 
/Cab (27) that emerges due to cross-correlations between 
two disorders a and b. 

Eq. (4.3) is universal and applicable for both homoge- 
neous RH-RH and mixed RH-LH stack-structures. The 
only difference between these cases is the sign in the un- 
perturbed phase shift fb = ±ujdb/c. This affects the 
value (3.8) of the Bloch phase 7 and changes the sign at 
the third cross-correlation term. 

For both homogeneous RH-RH and mixed RH-LH 
stack-structures, the Lyapunov exponent typically obeys 
the conventional frequency dependence, 

A = d/Li oc oc ui 2 when u> -> 0. (4.4) 

However, specific correlations in the disorders of the re- 
fractive indices taken into account in Eq. (4.3), may re- 
sult in a quite unusual w-dependence. In this respect, 
of special interest are long-range correlations leading to 
significant decrease or increase of the localization length 
Li oc (oj) in the predefined frequency window. Due to these 
correlations one can enhance or suppress the localiza- 
tion in the systems with compositional disorder (see, e.g., 
Refs. 9, 10). 

The expression for the Lyapunov exponent X(oj) man- 
ifests en emergence of the Fabry-Perot resonances asso- 
ciated with multiple reflections inside a or 6 slabs from 
the interfaces. These resonances occur when the wave 
frequency ui meets the conditions, 

lo/c = s a Tr/ d a or w/c = s b w/d b , s a ,6 = 1,2,3, 

(4.5) 

At the resonances the factor smip a or sin<y9 in Eq. (4.3) 
vanishes, resulting in the resonance increase of the lo- 
calization length Li oc . When only one type of the basic 
layers is disordered i.e., Eq. (4.3) contains only one corre- 
sponding term, the resonances give rise to the divergence 
of Li oc (lu). In the special case when the ratio between 
d a and db is a rational number, d a /db — s a /s , some 
resonances from different types of layers coincide. This 
also leads to the divergence of the localization length. 
The unexpected feature of these resonances is that they 
are quite broad due to vanishing of smooth trigonomet- 
ric functions. As was shown in Refs. 27, 28, 32, 33, in 
the case of positional disorder the terms entering the Lya- 
punov exponent and related to auto-correlations between 
the same type of slabs display the Fabry-Perot resonances 
associated with the other type of layers. On the contrary, 
in the case of compositional disorder the corresponding 
terms in Eq. (4.3) manifest the Fabry-Perot resonances 
belonging to the same type of slabs. 

As an example, let us consider the particular case of 
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FIG. 1: (color online) Lyapunov exponent versus frequency. 
Top: the mixed RH-LH material. Bottom: the RH-RH media. 
For both cases o a ~ cri ~ 0.006, d a — 0.6 , db = 0.4 , c = 1 and 
the length of sequences is N = 10 6 . Smooth curve corresponds 
to Eq.(4.7). 



the white-noise disorders for a and b slabs, 

a 2 ab = Q, K a {k)=K, b {k) = l. (4.6) 

In this case the Lyapunov exponent for both homoge- 
neous RH-RH and mixed RH-LH stack-structures takes 
a quite simple form, 

\= -^L = ha 2 a sin 2 va + v 2 b sin 2 <pb)- (4.7) 

Numerical results perfectly confirm this dependence, 
see Fig. 1. The more detailed comparison for iodj c <C 1 
and ujd/c 3> 1 also shows a nice correspondence. The 
data shown here are obtained with the use of Eqs. (3.11) 
without any approximation. One can see that for a long 
enough sample and weak disorder the analytical expres- 
sion (4.7) perfectly corresponds to the data, apart from 
random fluctuations. For each case only one realization 
of the disorder was used, and the fluctuations can be 
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smoothed out by an additional ensemble averaging. In 
order to see whether our predictions can be applied in 
experiment, we also used a quite short sample, N = 100, 
and strong disorder, a a as Ob « 0.3. The result shows 
that the analytical expression is also valid in average, 
however, the fluctuations around the smooth analytical 
curve are larger. 



V. MIXED RH-LH STACK WITH VANISHING 
BLOCH PHASE 

As we already noted, the expression (4.3) for the Lya- 
punov exponent is valid for both the homogeneous RH- 
RH and mixed RH-LH bi-layer array when the Bloch 
phase (3.8) is essentially different from zero. In this case 
the unperturbed map (3.10) does not degenerate, and 
nothing special arises for the evaluation of the Lyapunov 
exponent. Physically this is related to the fact that the 
unperturbed phase 9 n rotates when passing the sample of 
length N — L/d resulting in a homogeneous randomiza- 
tion. This happens everywhere provided the Bloch phase 
7 are irrational with respect to 2tt. Note that in this case 
the ^-distribution can be regarded as homogeneous even 
without the disorder. 

The situation is completely different in the case when 
the wave phase 7 = ip a + ipb vanishes after passing the 
unit (a, b) cell. This happens in the ideal mixed RH-LH 
stack with d a = db because after passing the RH a-layer, 
the phase shift is <p a = ud a /c, however, it is exactly can- 
celed by another shift, ipb = —Lud a /c = —ip a in the next 
LH fe-slab. As one can see, the circle (3.10) presenting 
the unperturbed map, degenerates into a single point in 
the phase space (Q,P), and, therefore, the unperturbed 
phase distribution is simply delta-function. This means 
that with a weak disorder, the phase distribution should 
not be expected as homogeneous one. 

In what follows, for simplicity we consider the ideal 
mixed stack whose refractive indices, n a and rib, ar e per- 
turbed by two independent white-noise disorders with the 
same strength [23, 29], 

o-l=o-l = a 2 , o- 2 ab = 0, K a {k)=K b {k) = l; (5.1a) 

d a = db, i.e., (p a = —<fb = and 7 = 0. (5.1b) 

The phase distribution p{9) can be found in the similar 
way as was described in Refs. 31, 36. The starting point 
is the exact recurrent relation (3.11b) between and 
n . By expanding this expression up to the second order 
in perturbation parameters rj a ^{n), one obtains 

0n+i-0n = -[Va(n)-Vb(n)}v(9 n )+a 2 v(0 n )v'{0 n ), (5.2) 
where we introduced the function 



v(9) = ip + sinc/3cos(20 — ip). 



(5.3) 



In deriving Eq. (5.2), we kept the linear terms propor- 
tional to T) a (n), r)b(ri) and substituted the terms r) 2 (n), 



rj 2 (n) by (rj 2 (n)) = (rj 2 (n)) = a 2 . Also, we explicitly took 
into account the condition (r] a (n)r)b(n)) = that directly 
follows from Eqs. (2.3) and (5.1a). This approximation 
is sufficient in order to obtain the distribution of phases 
9 n in the second order of perturbation. As a result, the 
expression (5.2) takes the form of the stochastic Ito equa- 
tion [38] which can be associated with the Fokker-Plank 
equation for the distribution function P(9,t) (see also 
Ref. [37]), 

The next step is to obtain the differential equation for 
the probability density p(6). In the case of weak disorder 
the corresponding Fokker-Plank equation for the distri- 
bution function P(6,t) has a relatively simple form, 



8 2 P f) 2 



2 d_ 

de 



P(6,t) 



dv 2 {9) 
dO 



(5.4) 

Here the "time" t is the same as the length A of a sample 
along which the evolution of phase 9 occurs. 

Since we are interested in the stationary solution of 
this equation, p{9) = P(6,t — > 00), the equation for p(9) 
reads 



d9 2 



[p(0)v 2 (0)] 



1 d 

2d9 



= 0. 



(5.5) 



Here the dependence on the variance a 2 has disappeared 
due to the rescaling of time, t — > a 2 t. Therefore, in this 
approximation the phase probability density p(9) does 
not depend on the disorder variance a 2 . Note that the 
only function v(9) entering the diffusive equation (5.5), 
is periodic with the period n. Consequently, its solution 
p{9) should be also periodic with the same period. In 
addition, p{9) should satisfy the normalization condition, 



p(0 + n)=p(0), 



^ d9p{9) 
Jo 



= 1. 



(5.6) 



One can easily obtain that the solution of Eqs. (5.5) and 
(5.6) is 



p(9) = J/v(0), J = - J ip 2 - si 

7T v 



sm 2 ip . 



(5.7) 



Our results indicate that the phase distribution p(0) 
strongly depends on the phase shift ip, and is highly in- 
homogeneous in the limit ip = ip a = u>d a /c <C 1, i.e., for 
small values of the wave frequency lo. Some examples 
of the distribution function p(0) for different values of 
ip are shown in Fig. 2. This figure clearly demonstrates 
that with a decrease of u> the distribution p(0) starts to 
be very sharp in the vicinity of = ir/2. 

It is worthwhile to note that the situation for the van- 
ishing value of 7 is somewhat similar to that emerging 
for the Anderson and Kronig-Penney models at the band 
edges. Indeed, in these models the unperturbed Bloch 
phase 7 also vanishes when approaching the band edges. 
This leads to a highly non-homogeneous distribution of 
perturbed phase, and, as a result, to a non-standard de- 
pendence of the Lyapunov exponent on the model pa- 
rameters. 
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FIG. 2: (color online) Stationary distribution p(8) for various 
values of uj = 95 in rescaled units d = 1, c = 1. Broken curves 
correspond to numerical data with an ensemble average for 
N = 10 6 , 10 7 , 10 8 for u> = 27r/5,27r/7,27r/15, respectively. 
Smooth curves present the analytical expression (5.7). 

In the considered model of the mixed RH-LH array, the 
crucial difference is that the effect related to the value 
7 = emerges independently of the frequency uj. This 
is in contrast with the case of Anderson and Kronig- 
Penney models for which the zero Bloch phase occurs 
at band edges only, therefore, for specific values of fre- 
quency. Thus, one can expect that for mixed RH-LH 
bi-layer stacks with the specific condition d a = db, the 
dependence of the Lyapunov exponent on the model pa- 
rameters has to be highly non-trivial. This fact can be 
seen from Fig. 3 which shows how the dependence of the 
Lyapunov exponent on frequency uj is affected by the re- 
lation between d a and db- 

The data in this figure clearly display that when d a 
approaches db, the standard dependence A oc uj 2 that is 
known to be a generic case for small uj, alternates by 
a very unusual dependence A oc uj 6 . The latter result 
was found numerically in Ref. 23 and later on, was dis- 
cussed in Refs. 29, 30. From Fig. 3 one can also see that 
for the homogeneous RH-RH stack-structure the conven- 
tional dependence (4.4) remains valid even in the specific 
case of d a = db. As was shown above, in this case the 
Lyapunov exponent A is described by the generic expres- 
sion (4.3) that for homogeneous RH-RH array is valid for 
any relation between d a and db- 

The Lyapunov exponent can be derived according to 
definition (3.12) and exact Hamiltonian map (3.11a) for 
the radius R„ . Its expansion within the second order of 
approximation in the disorder rj a ^{n) gives 

A = 2cr 2 sin </j(cos(20„ - <p)v{d n )) . (5.8) 

The averaging in this expression has to be performed with 
the distribution function p(9) determined by Eqs. (5.7) 
and (5.3). Since the denominator v(9) in Eq. (5.7) is 
the same as the coefficient in Eq. (5.8), we come to the 



remarkable result that the Lyapunov exponent (5.8) van- 
ishes for any value of the phase shift tp ! This means that 
in order to derive the non- vanishing Lyapunov exponent, 
one has to obtain the expressions for both the phase dis- 
tr 
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FIG. 3: (color online) Inverse Lyapunov exponent versus the 
rescaled frequency uj with the variance of disorder a 2 = 0.02 
and sequence length N = 10 6 . The data are obtained with an 
ensemble averaging over 100 realizations of disorder; a) mixed 
RH-LH array with d a = db, b) the same for d a = 0.99(4, c) 
the same for d a = O.ldb, d) homogeneous RH-RH array with 
d a = db- For comparison, the dot-dashed lines show two 
frequencies dependencies discussed in the text. 



of perturbation, by expanding them up to the fourth or- 
der in disorder which is not a simple task. Thus, one can 
expect that the Lyapunov exponent for the ideal mixed 
RH-LH stack with d a = db, should be proportional to c 4 
in contrast with the conventional quadratic dependence, 
A oc a 2 . Our extensive numerical and analytical studies 
confirm this expectation. 
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